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We study resolved sideband laser cooling of a one-dimensional optical lattice with one atom per 
site, and in particular the effect of the dipole interaction between radiating atoms. For simplicity, 
we consider the case where only a single cooling laser is applied. We derive a master equation, and 
solve it in the limit of a deep lattice and weak laser driving. We find that the dipole interaction 
significantly changes the final temperature of the particles, increasing it for some phonon wavevectors 
and decreasing it for others. The total phonon energy over all modes is typically higher than in the 
non-interacting case, but can be made lower by the right choice of parameters. 
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I. INTRODUCTION 

Over the recent decades, a variety of laser cooling tech- 
niques have been developed PHI] and implemented [5HD] 
for cooling trapped atoms or ions to very low tempera- 
tures. The resulting ultracold atoms can be made into 
a Bosc-Einstein condensate [SJ [D] . They can also be put 
into an optical lattice (an optical standing wave forming 
a periodic potential [TUlll2| ). where they can be used for 
quantum simulation of condensed matter systems [12], 
or potentially for quantum information processing [ID] . 
However, laser cooling does have limitations, which often 
require it to be followed by the loss of ~ 99 — 99.9% of 
the atoms by evaporative cooling [5] [6] . 

We here consider resolved sideband cooling [TJ 0], 
which is mostly used on ions [13] but has been applied 
to atoms in an optical lattice [21 [15] . This cooling tech- 
nique works on strongly confined particles: the trapping 
frequency v needs to be large compared to the sponta- 
neous decay rate T of the internally excited state. This 
allows the quantized motional states (with spacing hv) 
to be individually resolved. We can hence tune the driv- 
ing laser to drive mainly a transition to the next lower 
motional level (Fig. [I]), removing one motional quantum 
per cooling cycle. 

Equilibrium is reached when this cooling is balanced 
by heating. For a single particle, heating can occur ei- 
ther by off-resonant excitation or by the random recoil 
of spontaneous emission 1 . As confined particles can 
only be heated in whole quanta hv, heating is rare if 
hv is large compared to the widths of these heating pro- 
cesses, hT and E R = h 2 k 2 ,/2m respectively (here k is the 
transition wavenumber and m is the atomic mass). This 
gives an equilibrium temperature of ~ hT 2 /v, compared 
to ~ hT for free space Doppler cooling p]. Resolved 
sideband cooling also requires only one laser frequency, 
rather than the spread of frequencies required for narrow- 
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FIG. 1: Energy level diagram of our system. In the left site we 
show the conventional sideband cooling cycle: for red detun- 
ing A « v, the resonant excitation process is to the level with 
one fewer phonons (motional quanta), while in a deep lattice 
the dominant spontaneous decay process does not change the 
phonon number. On the right, we show an example of dipole- 
dipole interaction: the photon emitted by the middle atom 
is reabsorbed by the right-hand one, moving the excitation 
from the middle to the right-hand site. In this example the 
phonon states of both atoms are unchanged, but we will see 
later that this is not always the case. 



line (hT < En) Doppler cooling [2]. 

When laser cooling is applied to many closely spaced 
particles, the spontaneous decay rate is modified by in- 
terference between the emission from different particles, 
which is called superradiance [16-18 \. It is also possible 
for a photon emitted by one atom to be reabsorbed by 
another [19]. These two effects together are the dipolc- 
dipole interaction we consider. This interaction exists 
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FIG. 2: (Color online) Physical layout of our system. We 
consider a one-dimensional array of tightly confined atoms, 
one per lattice site. For simplicity we assume only one cool- 
ing laser, and hence can cool the atomic motion in any one 
dimension, either along the lattice (which we call x) or across 
it (y or z), but not all three at once. As in the single parti- 
cle case, cooling all three dimensions would require multiple 
lasers. 



between any particles emitting electric dipole radiation 
(the dipole is the matrix element of the transition), and 
becomes significant when the distance between them is of 
the order of the transition wavelength. It modifies the de- 
cay rate and propagates electronic excitations from atom 
to atom, and has no effect when no electronic excitations 
are present. It is hence qualitatively distinct from static 
interparticle forces, such as the Coulomb interaction be- 
tween ions [20] or the static-dipole interaction between 
aligned polar molecules [21] . which propagate motional 
excitations (phonons). Previous treatments of this in- 
teraction have considered either the superradiance part 
only [22] , the reabsorption part only [23] , or a long range 
approximation [TjJJ [23] which is valid for a moderate den- 
sity gas but breaks down at the half-wavelength spacing 
typical of an optical lattice. 

In this paper we consider a deep one-dimensional op- 
tical lattice with one atom per site, as shown in Fig. [2] 
and include the full dipole-dipole interaction. We solve 
for the steady state by transforming to momentum space 
and using a diagrammatic perturbative expansion in the 
Lamb-Dicke parameter rj = dE^Jv. We find that the 
dipole-dipole interaction may increase or decrease the 
steady state phonon number, depending on the laser pa- 
rameters. With optimized settings, the phonon number 
can be made ss 30% lower than would be possible without 
interaction. 



There are six sections in this paper. Section [IT] derives 
the relevant master equation for two-level atoms tightly 
confined by an optical lattice, including the dipole-dipole 
interaction. In Section |III| we simplify this equation for 
the case of an infinite one-dimensional optical lattice by 
transforming to momentum space. In Section [IV[ we in- 
troduce a Feynman-like diagrammatic perturbative ex- 
pansion, and use it to solve for the steady state phonon 
number distribution to lowest order in the Lamb-Dicke 
parameters rj a . The Appendix does the same expansion 
using non-diagrammatic methods. In Scction[Vj we eval- 
uate this steady state phonon distribution for particu- 
lar lattice setups, and find the parameter settings that 
minimize the phonon number. Finally, Section |VI| gives 
conclusions. 



II. THEORETICAL MODEL 

In this section, we derive the master equation of our 
system, including both the dipole-dipole interaction be- 
tween the atoms and the coupling between their elec- 
tronic and motional states. 

Our starting point is the dipole Hamiltonian, which 
includes energy terms for the electronic and motional 
(phonon) excitations of the atoms and for free photons, 
and interaction terms between the atomic dipole mo- 
ments and the electric fields of the cooling laser and the 
free photons. We then eliminate the photons to obtain 
the spontaneous decay and dipole-dipole interaction, us- 
ing similar methods to Refs. [TBHT5] . However, they con- 
sider atoms in fixed positions, while we allow the atoms 
to move within their lattice sites. We hence obtain addi- 
tional terms in the dipole-dipole interaction, describing 
processes that change the motional state as well as the 
electronic state. We derive explicit expressions for the 
corresponding coupling constants up to second order in 
the Lamb-Dicke parameters. 



A. The Schrodinger picture 

Our system consists of N two-level atoms in a deep 
ID optical lattice with one atom per site, as shown in 
Figs. |I|2[ We denote the two internal levels of atom i 
by |0) i and |l) f , and the transition energy between them 
by hujQ. We then define the atom raising and lowering 
operators 

S? = |1>«<0| and S* = |0>„ <1] . (1) 

We assume that the lattice is sufficiently deep that 
tunneling between neighboring lattice sites is negligible 
(i.e. atom i stays in site i), and each site can be ap- 
proximated by a harmonic well. We denote the trap- 
ping frequency of this harmonic well by v a , which we 
assume to be the same for both internal states, but allow 
to be anisotropic, where a = x, y, z denotes the direction. 
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However, unlike the usual Bose-Hubbard model [TU], we 
do not assume that the atoms are in the lowest motional 
band of the lattice, but allow them to be in higher mo- 
tional states. The excitations to these higher motional 
states are the phonons, with creation and annihilation 
operators B\ a , B ia for a motional excitation of atom i in 
the direction a. 

For simplicity wc consider the case of a single plane po- 
larized cooling laser, with peak electric field El, wavevec- 
tor kL, and frequency wl- To describe the spontaneous 
decay and dipole-dipole interaction, we also need to in- 
clude the free radiation field, with photon annihilation 
operators OkA labelled by a wavevector k and a polariza- 
tion index A = 1,2. 

We treat the atom-field interactions in the electric 
dipole approximation, i.e. as the dot product of the 
atomic dipole moment 

D, = Doi5 t +H.c. (2) 

and the electric field. The classical laser field is 

E L (x,t) = E L cos(k L • x - urf) , (3) 

while the free photons have electric field operator 



E(x) = ]Ti 



kA 



2e hL 3 



e k A e Ikx a kA + H.c. , (4) 



where L 3 is the quantization volume and e^A the polar- 
ization vector. 

Putting all these together, the Hamiltonian of the sys- 
tem in the Schrodinger picture is 



JV 



/v 



H = J2 huj o S l S i+J2 e ' Di - [E L (x^i)+E( Xi )] 



i=i 



i=i 



N 



^ fi^fe « k A a kA + ^ Hu a B\ a B ia , (5) 

i—l ct=x,y,z 



kA 



where Xj is the position of atom i. 



B. The interaction picture 

We next eliminate the time dependence of the laser 
term, by moving to the interaction picture with respect 
to 



N 

H = ^2 SjS l + ^2 a kA akA , (6) 

»=1 kA 

and applying the rotating wave approximation with re- 
spect to wq (i-e. neglecting terms e 1 "' with oj > uiq). This 



gives the interaction picture Hamiltonian 

N N 



i=l L 



i=l 



kA 



U>k T\ ,* .-ik-x; i(ui t -u L )(g t 



N 



H.c. +Y, hv a B\ a B ia , 



(7) 



2—1 ct— x,y,z 



where we define the laser detuning A = luq — cjl (i.e. red 
detuning is positive) and the Rabi frequency 



n = ^ D i ■ E L . 



(8) 



We choose the relative phase of |0), |1) so as to make f2 
real. 



C. The master equation 

As our system is not in an optical cavity, the free pho- 
tons a k A quickly leave the system, and are absorbed by 
the walls of the apparatus on a time scale At ~ L/c, mak- 
ing our system an open quantum system. We hence need 
to use a master equation rather than a Hamiltonian to 
describe its evolution. For realistic sizes, ujq 3> 1/At^> 
v, A,f2,r (where V is the decay rate derived below), 
which allows us to use the method of Refs. [TH1 El H5J to 
derive a Markovian (memoryless) master equation. 

This method approximates the photon loss process by 
alternating free evolution under H\ for time At with the 
instantaneous loss of all photons, 



p(t + At) = Tr photons [Ui(t + At, t ) 

x{ P (t o )®|0) (0| pht }E# o +AMo) 



(9) 



where pit) is the (electronic and motional) density matrix 
of the atoms at time t, |0) ht is the vacuum state of the 
photon field, and Ui(to + At, to) is the unitary evolution 
operator under Hi from time to to to + At. 

Since At is short compared to the dynamical 
timescales of Hi, we can approximate this by expand- 
ing the time ordered exponential Uiito + At, to) — 

( _ to+At ~l 
Texp < — ^ J dtHiit) > up to second order in At: 



Ap 
At 



Tr 



photons 



to+At t 



to+At 




dtHi(t) 



(10) 



l 

hAt 



dt / dt'Hi{t)Hi(t') \ p(to) ® |0> (0| pht + H. 



to *0 
t +At to+At 



h 2 At 



dt / dt'H I (t)p(t o )®\0)(0\ pU H l (t') 



to to 
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where Ap = p(t + At)-p(t Q ). In the limit u > 1/At> 
v, A, fi, r, the above integrals can be evaluated and give 
a time-independent master equation, which can be con- 
veniently written as [TH1 Q~7] 



UH condP -pHl ond ] +K(p) 



(11) 



where i? C ond comes from the terms of ( 10 1 in curly brack- 
ets and 1Z(p) comes from the last term. 



D. The no-photon time evolution 



The term i? con d m Eq. (11) describes processes where 



no photons are lost from the system, cither because none 
are emitted, or because one is emitted but is reabsorbed 
before time At [El Hi] . As it appears in the master equa- 
tion (11 ) in a similar form to the — h [H, p] of a Hamilto- 



nian, but describes only evolution under the condition of 
no photon loss, this term is called the conditional Hamil- 
tonian. Unlike a normal Hamiltonian, it is not Hcrmi- 
tian, but norm decreasing; this represents the decreasing 
probability that this condition will continue to hold, and 
in the master equation is balanced by TZ(p). 

Evaluating (10) using the methods of OUT], we find 



H, 



cond 



N 

E 



^ e - ikLXi 5i + H.c. 



^^(x.-x,)^ 



3& 



(12) 



The terms of (12) describe respectively laser driving, 



detuning and spontaneous decay, phonon energy, and 
dipole-dipole interaction. All except the last are iden- 
tical to standard laser cooling. 

The last term contains the dipole coupling operator 



C(r) = -Te ik ° r 



1 

ik r 



1 



;i-i£>oi-fi 2 ) 



(k Q r) 2 i(fc r) s 



(l-3|D 01 -f| 2 ) 



,(13) 



where r = ||r||, f = r/r, fc = uj /c and D i = 
Doi/||Doi||. The imaginary part of C(xj — Xj) describes 
photons being emitted by atom i and reabsorbed by atom 
j, while the real part describes modifications to the de- 
cay rate due to interference between their emitted waves. 
Eq. ( 13 ) is exactly the same expression as the dipole cou- 

but they fix the atom 



pling constants in Rcfs. p!6H18] 
positions so their r is a classical vector, while we allow 
motion within each lattice site and our r is hence an op- 
erator. 



E. Lamb-Dicke expansion 

Eq. |l2"| ) is difficult to use as it stands because it con- 
tains transcendental functions of the operator Xj, namely 
exp(— ikL • Xj) and C(xj — x,). In this section we intro- 
duce the Lamb-Dicke expansion about the site centers, 
which allows these to be handled perturbatively in the 
deep lattice limit. 

As already mentioned above, we approximate each 
lattice site by a harmonic well with phonon operators 
B ia , B\ a . In terms of these operators, the position oper- 
ator X* is 



(14) 



where is the equilibrium position (site center), a — 
x, y, z labels vector components, and the Lamb-Dicke pa- 
rameter r\ a is defined by 



■q a = Jhkl/2mv a 



(15) 



where m is the mass of the atom. In the case of a deep 
lattice rj a is small, and we can hence use Eq. ( 14 ) to 



approximate exp(— ikL 1 x i) an d C(xj — x,) in Eq. (12) 
by Taylor expansion. Up to second order in rj a , we find 



e — ik L -Xj _ e -ik L Xi 



1-i h a r)a(Bia + B] a ) 



a—x,y,z 



E ^hJ^{B*+Bl){Bv + B, 



ill) 



a,0—x,y,z 



(16) 



where kL = Icl/IIcl | and k^ a is the component of this 
unit vector in the a direction, and 

C(Xj — Xj) = Cij 

+ E D U B ^+B] a -B la -B\ a ) 

a=x,y,z 

+ E E rf(B ja +Bl-B ia -Bl) 



a : p—x,y,z 



y-{B jP +B} p -R 



i0 



B 



(17) 



a/3 



Here the dipole coupling constants Cy, Df-, and E 
are obtained by differentiating C(r) (and absorbing the 
Va/ko factors), 



D a (r) 



tfc dC(r) 
k dr a 

32/ 



y ' " 2fc2 dr a drp 



(18) 



where r a are the components of r, then substituting the 
equilibrium positions Xj for the atoms, 



Cij = C(Xj — Xj) 



D'; 



E 



a/3 _ 



D a (Xj-Xi), 



(19) 
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Because they contain only equilibrium positions and 
not position operators Xj, these are numbers rather than 
operators. Finally we remark that C'ij — Cji, and hence 



D 



a/3 



„ l>j, ^d E^f = Ej. In Fig 
coupling constants against the distance 
the atoms, for a few relevant angles between Xj 
the dipole moment Dqi- 



3 we plot these 
Lj— Xj| between 
X, and 



Substituting Eqs. (16) and (17) back into the condi- 
tional Hamiltonian (12), we see that in zeroth order in 
77 a , the dipole-dipole coupling constants Cij are the same 
as for atoms in fixed positions [T6HT8] , while the motional 
states are unaffected. The first and higher order terms 
couple the electronic and motional states of the atoms: 
when atoms exchange a photon, they may also create or 
destroy one or more phonons, with coupling constants 
Dg and Eg. 



F. The effect of spontaneous emission 



The term TZ(p) in Eq. (11) describes processes where 
photon loss does occur; it is called the reset operator as 
it "resets" the system state after such a loss [35] • To 
evaluate this expression we proceed as in Refs. [TrJl IT?] , 
and again obtain something similar to the fixed-positions 
expression ReC(x :) — x^SipSj with the positions be- 
coming operators. However, here ReC(xj — x,) does 
not simply act from the left, but becomes a superop- 
erator (i.e. contains operators acting from both sides). 
In Lamb-Dicke expansion, each phonon operator acts 
from the same side of p as the same atom's electronic 
(de)excitation operator, i.e. Bi a + B\ a on the left and 
on the right. Up to second order in rj a , this 

gives 

N 

i=l j^i 

+ E E ReD%[S iP S}(B ja + Bl) 

j^i a.—x,y,z 

-{B ia + B\ a )S iP S]] 

N 

+ E E Rc E t [ {Bia + Bl) (B ip + Bl)Si p S] 

i,j=l ot,/3=x,y,z 



-2(B ia + Bl a )S iP S}(B jlB + B^ 



S iP S}{B ja + B} a )(B j0 + B} )] 



(20) 



As in the case of fixed positions, the spontaneous emis- 
sion rate can be higher (superradiance) or lower (sub- 
radiance) than its non-interacting value of T per ex- 
cited atom, depending on the distance between the 
atoms. With motion included, there is also the possi- 
bility that a spontaneous emission creates or destroys 
phonons, at rates ReDfj, ReE"?. The single particle 
recoil heating is included as the i = j term ReE^f = 
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FIG. 3: (Color online) Dipole coupling constants C(X), 
D a (X.), and E aa (X), as a function of distance |X|, for along- 
lattice motion (black) and across-lattice motion (blue). Solid 
lines are $ = (Doi ■ X = sini? for along-lattice motion or 
cosi? for across-lattice motion, cf. Fig. [2j). Dotted lines are 
i? = 0.56 (along-lattice) or 0.35 (across-lattice), which we will 
see later are the optimum values at lattice spacing a = 0.5 A. 
For D and E the component plotted is the relevant one for 
cooling that motion direction: a = x for along-lattice motion 
or a = y for across-lattice motion. The vertical lines mark the 
lattice positions for a = 0.5A (blue dotted) and a = 0.744A 
(black solid). 



ion 



'01/3 



25ap ) i where 5 a p is the Kro- 



necker delta. 

In the limit of infinitely stron g co n finem ent (rj a — > 0, 
and hence Df 



E 



a/3 



■ij, ~ij _ ' 0), Eqs 
known master equation for fixed c 



12), (20) reduce to the 
lpole interacting atoms 



G 



III. FOURIER MODE REPRESENTATION 

The lattice has discrete translational symmetry, and 
we can hence simplify the above Hamiltonian ( 12 1 and re- 



set operator (20) by transforming to Fourier space, which 



we do in this section. 

A. Mode operators 

The Fourier mode operators s p , b pa are labelled by 
a wavevector p, which runs over the first Brillouin zone 
—n/a < p a < 7r/a, where a is the lattice spacing. They 
are defined in terms of the single site electronic excitation 
and phonon operators Si, Bi a by 

1 N 

v i—1 

i N 

b pa = 77^I>~ ip ' XlB »' ( 21 ) 



where N is the total number of atoms and we assume 



periodic boundary conditions. The creation operators s 



b pa are given by the Hermitian conjugates of Eq. ( 21 ) . As 
before, a denotes the direction of vibration of the atoms 
and can assume the values x, y, and z. 

These operators create or destroy collective excitations 



spread over all the atoms. The internal atomic excita- 
tions are hardcore bosons (no more than one can exist 
on any one atom), and hence s p , s p do not exactly obey 
normal bosonic commutation relations, but in the limit 
of small excitation density the difference can be neglected 
[26l 127] . (We cannot use the Jordan- Wigner transforma- 
tion to turn them into fermions, as their number is not 
conserved.) The collective phonons described by b pa are 
ordinary bosons, as in the deep lattice limit one atom 
can have any number of phonons. As we shall see below, 
the electronic excitations propagate by photon exchange, 
while the phonons do not propagate to zeroth order in 



B. The conditional Hamiltonian 



The definitions (21) can be inverted to give the sin- 
gle site operators Si, Bi a in terms of the Fourier mode 
operators s p , b pa . We now substitute this into the con 



ditional Hamiltonian Eq. (12), and make a Lamb-Dicke 



expansion up to second order in rj a , using Eqs. ( 16 )-( 18 1 
Introducing the notation 



r(o) 

co nd 



f (i) 

cond 



H 



(2) 



ad > 



(22) 



where the superscripts indicate the order of each term 
with respect to r\ a , we obtain 



H 



(0) 
cond 



H, 



(i) 

cond 



H, 



(2) 
cond 



hVNn(s kh + 4 J - Y \ hc{ v) 4> s p + E hiy c* b U b P* 



p . a 



p,Q p,q,a V 

- E 



rj a r] P nk La k L/3 (b pa + feL pQ )(6 q/ 3 + 6i q/ 3)(4 L + P + q + s k L - P -q) 



p,q>«,/3 

E ■ 

p.q,r,a,/3 



4VA 



"2N 



[e af) (r + q) - 2e afj {r) + e afJ {r - p)](b pa + 6L pa )(6 q/3 + ^ q ^)4+q s >--P 



(23) 



where addition of wavevectors is defined mod the Bril- 
louin zone. Here the coefficients c(p), d a (p), e a p(p) are 
given by 



c(p) = ^e- i P' x -C(X l ), 

i 

d a (p) = ^e-^LPpLi), 

i 

e a p(p) = Y c ~ iP ' X * EaP ^)- ( 24 ) 



These are the discrete Fourier transforms of the constants 



dj, Df n and in Eqs. rtlQh and <\18\, with the single 



atom terms h(A — ±T) absorbed into them by allowing 
i = j and defining 



C(0) = r + 2iA and D a (O) = 0. 



(25) 



They are functions of wavevector p, and also depend 
on the lattice spacing and orientation; they are plotted 
against p in Fig. (9Td) for a few example settings. 

The main advantage of changing into the Fourier space 
is that the zeroth order term (in rj a ) C(X) in the dipole- 
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dipole interaction becomes the dispersion c(p) of collec- 
tive electronic excitations, which allows this term to be 
treated non-perturbatively. The first and second order 
terms d a (p), e a p(p), which we do still have to treat per- 
turbatively, describe interaction processes where an elec- 
tronic excitation emits or absorbs one or two phonons. 
Because of the discrete translational symmetry, wavevec- 
tor (quasimomentum) is conserved mod the Brillouin 
zone, with the laser contributing its wavevector Icl to 
each excitation it creates. 

In the large system limit N — > oo, these Fourier sums 
converge in one dimension (except for a logarithmic di- 
vergence at the points p = ±ko) but diverge in two or 
more dimensions. While they can still be calculated for 
finite sizes, in the presence of such strong finite size ef- 
fects we cannot neglect the difference between periodic 
and open boundary conditions. Physically, this is be- 
cause in 2D and 3D photons have a significant probabil- 
ity of being reabsorbed after travelling a large distance; 



in an infinite 3D system there would be no outside for 
photons to escape to, and they would instead mix with 
the atomic excitations to form two polariton modes [25] . 
(In the infinite system limit there are also retardation 
effects neglected by our master equation, but for prac- 
tical experimental sizes the boundary is reached before 
these become significant.) In the following, we therefore 
restrict ourselves to one-dimensional lattices. 



C. The reset operator 



Similarly, the reset operator IZ(p) in Eq. (20) can now 
be expanded as 



K(p) = K i0) (p)+K {1) (p)+K (2 \p), (26) 



with 



TZ^(p) = ^Rec(p) Sp p4, 
p 

K {1) {p) = -nr T Imd "(q) s qP4+q( 6 P Q+6 -P«) _Imd «(P + ^)( 6 P« +6 -pa) s q^ 



p,q," 



K- {2) {p) = ft R.ee Q(3 (r + q)(6 pQ + bl pa )(b clfJ + fc^K-pPsJ+q 



p,q,r,a,/3 



2Ree Q/3 (r)(& pQ + &_ p Jsr- P K+ q ( 6 q£ + b -qfs) + R-ee a/3 (r - p)s r _ P K+ q ( 6 P« + bl pa )(b^ + b[ 



q/3' 



(27) 



Because e ipx + e tpx = 2 cos pec is real and e tpx — 
e -ipx _ 2isinpx is pure imaginary, even functions such 



as CijjE^f transform real part to real part (ReC., 



i j 



Rec(p)) while odd functions such as ReDg transform 
real part to imaginary part (Re-D™ — > iImd Q (p)). The 
zeroth order term describes the spontaneous decay of the 
electronic excitations, with a wavevector-dependent rate 
because of dipole-dipole interaction. The first and sec- 
ond order terms are decays that also create or destroy 
phonons. 



IV. STEADY STATE DENSITY MATRIX 



Under the master equation (11), the density matrix 



converges to a steady state, and for cooling we would like 
it to converge within a reasonable time to a state with 
the lowest possible temperature, or equivalently phonon 
number. In this section, we calculate this steady state 
using an open system version of Feynman diagrams. 



A. Open system Feynman diagrams 

In this subsection, we describe a modified form of Feyn- 
man diagrams that can be used on open systems de- 
scribed by a time independent Markovian master equa- 
tion. Our diagrams are similar to those of Keldysh 
[291 130] . but our starting point is a master equation rather 
than a system plus bath Hamiltonian. 

Our diagrams divide the Liouvillean superoperator 
C(p) — p into a propagation part and a vertex part, 
and treat the vertex part as a perturbation, in the same 
way as conventional Feynman diagrams divide the Hamil- 
tonian operator |31j . Propagation terms are those in 
-ffcond with a matching creation and annihilation oper- 
ator, i.e. — ^hc(p) s p s p for an electronic excitation and 
hv a b pa b pa for a phonon. (These lines represent propaga- 
tion through time as well as space, so it makes sense for 
a non-propagating particle such as our phonon to have a 
propagation line.) The remaining terms in i? CO nd and all 
terms in 7Z(p) are vertices. Fig. |4|a) shows the vertices 
appearing in our master equation. 

Propagation lines are labelled with whatever indices 
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FIG. 4: (Color online) Feynman diagrams. Solid orange line: 
electronic excitation, dotted green line: phonon, blue circle: 
laser (de)excitation, blue star: spontaneous decay. The black 
vertical line is used to keep track of which operators act from 
which side of p. (a) Table of allowed vertices, by number of 
phonon lines (order in 7/ a ) n and type of interaction. Vertices 
exist with arbitrarily high n, but for our purposes we need 
only consider up to 2nd order, i.e. n < 2. In the case of de- 
cay, the operator ordering shown is not the only one allowed: 
there is always one atomic excitation on each side of p, but 
the phonon(s) can be on either or both sides. Where no ar- 
rowhead is shown either propagation direction is allowed, (b) 
Example of a rate term (diagram symmetrical about the p 
line), the ordinary sideband cooling process, (c) Example of 
an interference term (diagram asymmetrical). Lines labelled 
with frequency E/h and wavevector p. 



are necessary to fully identify an operator, here p for an 
electronic excitation and pa for a phonon, and also with 
an energy E, which is generally not on that particle's 
dispersion. Each operator in a vertex must be connected 
to a matching propagation line, which may go to another 
vertex or leave the diagram. The energy E is conserved 
at vertices (total in=total out), and a particular master 
equation may have additional conserved quantities, e.g. 
wavevector if it is translationally invariant (such as ours) . 
Only particles appearing as operators in the master equa- 
tion appear as lines in the diagram; particles eliminated 
during the master equation derivation (here, decay pho- 
tons) do not. Classical entities such as the laser field also 
do not appear as lines, but may contribute momentum, 



here Icl at each laser vertex. 

A master equation contains operators acting on both 
sides of the density matrix p, and we hence need to keep 
track of which operators act on the left of p and which on 
the right, which we do with the black vertical line in the 
diagrams. Propagation lines and vertices from H con dP 
are on the left of this p line, while those from pH^ cond 
are on the right; a line or vertex on the right hence has 
the complex conjugate coefficient of an identical one on 
the left. Reset operator vertices contain operators acting 
on both sides of p, and hence appear only on the p line, 
with propagation lines attached to them from both sides. 
Propagation lines cannot cross the p line without going 
through a reset vertex. 

We choose to make the arrow on a propagation line 
point forward in time, i.e. from the vertex where that 
excitation is created to that where it is destroyed. Since 
s_ |0) = |l p ) but (l p | s p = (0|, this is the usual creation 
operator (of the vertex) to annihilation operator on the 
left side of the diagram, but annihilation operator to cre- 
ation operator on the right side. With this convention, 
energy/momentum conservation at a reset vertex means 
equal energy /momentum on each side of p. (The alterna- 
tive of creation operator to annihilation operator on both 
sides would also be self-consistent, and would make the 
conservation laws in=out everywhere, but would make 
the arrows point backwards in time on the right.) 

To evaluate such a diagram, we replace each vertex 
by its coefficient (in the master equation, i.e. including 
the ±i/h where applicable) and each internal line by its 
propagator l/(— coefficient =F where the — applies 

to the left side of p and the + to the right. Lines en- 
tering or leaving the diagram have on-resonance energy. 
As in conventional Feynman diagrams, for tree diagrams 
the external lines' energies (and momenta, if conserved) 
determine all the internal ones by conservation, while for 
diagrams containing loops it is necessary to integrate over 
the undetermined energies and momenta. For example, 
ordinary sideband cooling is represented by the diagram 
Fig. Kb), which has the value 



Rec(p + k L ) 




|c(p + k L )-2n/ a | 2 



(28) 



where p is the phonon momentum (which is the same on 
both sides by conservation). 

The resulting value is that diagram's contribution to 
the transition rate (not amplitude) between its start and 
end states, which are given by the particles entering and 
leaving it. Symmetrical diagrams such as Fig. 4](b) give 
the rate of the process appearing on both sides, while 
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asymmetrical diagrams such as Fig. [4^c) give the inter- 
ference term between their left and right processes. In 
the case of coherent processes (i.e. not involving any re- 
set terms), the left and right sides are disconnected, and 
the transition rate is simply the product of the left and 
right amplitudes. As already explained, amplitudes on 
the right are complex conjugated, so this is consistent 
with rate = | amplitude! 2 . 



B. Heating and cooling rates 



V(p> = v-- y = 



In this section we use these diagrams to calculate the 
heating and cooling rates. 

The lowest order (in the Lamb-Dicke parameters rj a ) 
processes that change the phonon number increase or de- 
crease it by 1 on both sides of p, as changing it on only one 
side is forbidden by energy conservation. This involves 
2 external phonon lines (one on each side) which must 
both end on vertices, and is hence at least 2nd order in 
7] a . In each case there are 6 such diagrams which must be 
summed to obtain the rate, shown in Fig. [5] for creating 
a phonon (heating); for destroying a phonon (cooling), 
they are the same except with all phonon lines reversed. 
For simplicity we restrict to the fully anisotropic case 
(three different phonon frequencies v x , v y , v z ), where en- 
ergy conservation requires that a be the same on both 
sides. 

Evaluating these diagrams gives 



-2Q 2 Rc e QQ (k L =FP) 



|c(k L 



i] a nk La ) Rec(k L =FP) 



|c(k L Tp)±2u/ Q | 2 
f2 2 |d Q (k L ) - d Q (k L t p)| 2 Rec(k L t P) 



-2Re 



|c(k L )| 2 |c(k LT p)±2i^| 2 
in 2 [d a (k L ) - d a (k L =p p)]Imrf a (k L =p p) 
|c(k L )| 2 [c(k L Tp)±2i^] 



i» 2 ?7 a fc La [d a (k L ) -d a (k L T-p)]Rec(k L j-p) 
c(k L )|c(k L Tp)±2i^| 2 

^ 2 ??afcLaImrf a (k L =F P) 

' c*(k L )[c(k L TP)±2ii/a] 



where p is the phonon momentum, which must be the 
same on both sides of p by conservation, and the up- 
per sign applies to phonon creation (i.e. heating) and 
the lower sign to phonon destruction (cooling). The 2Re 
on the 3 asymmetric diagrams is because we also need 
to include their mirror images, which have the complex 
conjugate value. 

At the same order in the r\ a , there are also a self-energy 
correction S Q (p) to the phonon line given by diagrams 
with one ingoing and one outgoing phonon line on the 
same side of p, and a vacuum term V given by diagrams 



0(r, 4 ) 



FIG. 5: (Color online) Feynman diagrams for the heating rate 
Rha (p) , as evaluated in Eq. 29 The cooling rate diagrams are 



the same except that all phonon arrows are reversed. As in 
Fig. [4] solid orange lines are atomic excitations, dotted green 
lines are phonons, blue circles are laser driving, blue stars are 
spontaneous decay. The black dotted circle denotes any al- 
lowed diagram with the given external lines. The line under 
each diagram is the color used to plot its value in Fig. [9^b) . 
The symmetrical diagrams (first 3) are rate terms; the asym- 
metrical diagrams (last 3) are interference terms, and are to 
be read as including both mirror-image versions. 



with no external lines. These are evaluated in the Ap- 
pendix (Eqs. (A.34), (A.35)). Together, these give an 
effective master equation for the phonons, 



P 



Vp 



p . a 



[(S a (p) - iv a )b pa b va p 



+(E* (p) + w a ) P bl a b va + R ha (p) 
+ R ca {p)b pa pbl a ] . 



DQ P ^PO 



(30) 



This conserves probability as Re S Q (p) = — \ (Rha (p) + 
-Rcq(p) and V = -Epa^ho- It also agrees with the 



(29) steady state equation (A. 33 1 calculated by direct pertur 



bative expansion in the Appendix. 



C. Steady state phonon distribution 

We now use these equations to calculate the steady 
state of the phonons. We also introduce the mean phonon 
number as a measure of heat content, as the steady state 
does not have a definite temperature. 



Eq. (30) does not couple different phonon modes pa, 



so the phonon steady state p^ h \ is a product state over 
modes. As ReSo,(p) < 0, phonon coherences decay, so 
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the steady state is a mixture (not a superposition) over 
phonon numbers, 



PpL = I1ES« \m pa )(m pa \ 



(31) 



Here the coefficient c„ lpQ is the probability of having m pa 
phonons in the mode pa. As the state is a product state, 
the phonon distributions of different modes are indepen- 
dent. 



Substituting ( 31 1 into ( 30 1 and setting p = 0, we obtain 



= m pa R ha (p)c r , 



[(* 



l)-Rha(p) + m pa R CQ (p)} c„ 



+ (m pa + 1) R ca (p) c„ lpQ+ i 
One solution of this equation is 

Rh a (p) 
R ca {p) 



(32) 



(33) 



and it is unique up to normalization because once cq is 
chosen the others can be determined iteratively. 

As expected for a steady state, heating and cooling 
are in detailed balance: the transition rate from to to 
m + 1 phonons is equal to the transition rate from m + 1 
to to. A normalizable steady state only exists for the 
modes where the heating rate Rh a (p) is smaller than 
the corresponding cooling rate R ca (p). If not, we obtain 
unbounded heating. 

We finally consider how to measure the heat content 
of this steady state, as in the next section we will want 
to minimise this. For each individual phonon mode pa, 
the state ( 33 ) is a thermal state with temperature 



T — 

± pa 



k B log [R ca (p)/R ha (p)] 



(34) 



However, the phonon state as a whole is not thermal 
because different modes pa have different temperatures. 
We therefore instead use the mean phonon number as our 
measure of heat content, as this is defined for any finite 
energy state, and when multiplied by hv a gives the mean 
phonon energy. For a single mode it is 



Rha{p) 



R ca (p) - Rha{p) 



(35) 



For multiple modes it can then be obtained by simple 
averaging; in particular, the mean phonon number per 
atom over all wavevectors is 



(TO) 



(36) 



We do not average over motion directions a, as in a single 
laser setup we expect cooling in one dimension (along the 
laser) and heating in the other two. 



V. EXAMPLES AND PARAMETER 
OPTIMIZATION 

In this section we evaluate the steady state phonon 
number for various parameter values, and in particular 
consider what settings minimize it. We find that the 
dipole-dipole interaction has a significant effect on the 
steady state phonon number, at lattice spacing a = A/2 
typically ~ 100% at the most strongly affected phonon 
wavevectors and ~ 20% averaged over all wavevectors. 
This effect can have either sign depending on the param- 
eters; in particular, there are settings that give a lower 
phonon number (m) ss than any non- interacting settings 
with the same lattice strength v. 



A. Ordinary (non-interacting) sideband cooling 

As a check on our calculation we first consider the case 
of atoms far enough apart to neglect their dipole-dipole 
interaction. In this case, we should obtain ordinary side- 
band cooling of each atom independently. No interaction 
means 



C i 



0, 



D a 



o, e: 



a/3 



Vi^j 



(37) 



The 



This leaves only the single atom terms Cu (Eq. (25)) 
and E°f, as = by definition (cf. Eq. Q). 
Fourier transforms are hence 

c(p)= C u = T + 21A, d a (p)=0, 

7 aP _ VaVP T 



e Q(3 (p) = 



10 



(d i q D 01/3 - 2<5 Q/3 ) ,(38) 



independent of p, where Doia is the component of the 
dipole unit vector Dqi along a. 

In ordinary sideband cooling, the cooled motion a is 
parallel to the laser wave vector (d = 0) and hence 
k^ a = 1 and Doi Q = for this motion. Substituting 



into Eq. ( 29 ) , we obtain the heating and cooling rates 



R ha = vl^r 



R ca = 7f a n 2 T 



1 



+ 



r 2 + 4(A + v a ) 2 5(r 2 +4A 2 ) 
1 2 



r 2 +4(A-!/ Q ) 2 5(r 2 +4A 2 ) 



(39) 



independent of p. These are the standard results for 
single particle sideband cooling [13l [32] . 

In this case, only the first 2 diagrams in Fig. [5] are 
nonzero. The first term in each of these equations rep- 
resents a phonon being created (heating) or destroyed 
(cooling) when the laser excites an atom, followed by 
spontaneous decay with no change in phonon state. Cool- 
ing is resonantly enhanced when A w expected 
from Fig. [T] The second term represents laser excitation 
with no change in phonon state, followed by spontaneous 
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FIG. 6: (Color online) Contour plots of the steady state 
phonon number {m) sa , averaged over wavevectors, as a func- 
tion of the laser detuning A and the angle $ between the 
atomic motion and the laser direction, (a) Non-interacting 
case, (b) Across-lattice motion, a = 0.5A. (c) Along-lattice 
motion, a = 0.5A. (d) Along-lattice motion, a = 0.744A. In all 
cases v — 10F and the atomic dipole moment Doi is perpen- 
dicular to the wavevector k.L of the incoming laser. All cases 
are symmetric about $ = 0; to save space, only positive # are 
shown here. The green shaded areas are where the phonon 
number is lower than its minimum non-interacting value. The 
blue dots are the minima, and the blue lines the sections plot- 
ted in Fig. [7] The contour interval is 10~ 4 phonons per atom. 



decay with phonon creation or destruction. This is re- 
coil heating; as it has equal heating and cooling rates, it 
alone would have an infinite equilibrium temperature, so 
it has a heating effect on any finite temperature state. 

Fig. [6^ a) is a contour plot of the steady state phonon 
number (m) ss as a function of the laser detuning A and 
the angle d between the cooling laser and the motion 
(cf. Fig. [2| . As expected, the lowest steady state phonon 
number is obtained for i3 = and A close (but not exactly 
equal) to the phonon frequency v a . 

Since we consider a one laser setup, the atomic motion 
/3 orthogonal to the laser wave vector is not cooled, but 
is still subject to recoil heating. For this motion we have 
= 0, and hence 



R 



h/3 



R, 



C0 



5(r2+4A2) \ 2 ~ D H 



(40) 



Both the cooling and the heating rate are the same for 
this motion. This means the mean phonon number (m) ss 
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FIG. 7: (Color online) Steady state phonon number (m) ss , 
averaged over wavevectors, against i9. The detuning A is in 
each case chosen so the plot passes through the minimum 
of phonon number (forming the sections of Fig. [6] along the 
blue lines). Dotted green: no interaction (a = 00). Solid 
blue: across-lattice motion, a — 0.5A. Solid black: along- 
lattice motion, a — 0.5A. Dotted black: along-lattice motion, 
a = 0.744A. 



for this motion diverges, i.e. unbounded heating, as ex- 
pected. A diagonally oriented laser could cool all three 
dimensions, but inefficiently as it can only be in reso- 
nance with one of them. Alternatively, three lasers could 
be used, one along each dimension. 



B. Interacting case: results 

We now move to the case with dipole-dipole interac- 
tion, and consider varying the lattice spacing a and the 
laser detuning A and orientation d, for both the along- 
and across-lattice modes (as defined in Fig. [2]). We re- 
strict to the case where the laser wavevector Icl is orthog- 
onal to the atomic dipole moment Dqi, which requires 
the least laser power for a given 17, and both are in the 
x-y plane. 

Figs. |6jb,c,d) are contour plots of the steady state 
phonon number (m) ss against the laser detuning A and 
angle for three choices of lattice spacing and orienta- 
tion. Fig.[7]plots the A = constant cross-section through 
the minimum of (m) ss , for the same cases. Fig. [8] plots 
steady state phonon number against lattice spacing. 

The numerical results presented here were all obtained 
at N = 200 with the dipole interaction cut off at that 
range (i.e. each pair of particles have one interaction 
term, not an infinite series of image interaction terms) 
and v = 10r. We have tried other N, v and found the 
same qualitative results for N > 5 and v > 5T (along- 
lattice motion) or > 2r (across-lattice motion). Increas- 
ing v decreases the steady state phonon number, for large 
v /r as approximately (m) ss ~ (r/z/) 2 , as in the non- 
interacting case. 
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FIG. 8: (Color online) Steady state phonon number (m) ss , av- 
eraged over wavevectors, against lattice spacing a, for along- 
lattice motion (black) and across-lattice motion (blue). Thin 
lines are standard sideband cooling A = v, # = 0, thick lines 
are A, i? chosen to minimize (m) ss . The dotted horizontal line 
is the non-interacting minimum. 



We first consider across-lattice motion with lattice 
spacing a — 0.5A (Fig. [6jb) ) . If we leave the settings 
at the non-interacting optimum of d — 0, A = 10.02r, 
the interaction increases the steady state phonon num- 
ber from 1.62 x 10~ 3 to 2.03 x 10~ 3 per atom. However, 
if we re-optimize A for the interacting system, we can 
reach a minimum of 1.51 x 1CP 3 phonons per atom at 
i9 = ±0.35, A = 10.12r. This is lower than the non- 
interacting minimum, i.e. dipole interaction can improve 
cooling effectiveness. In this case, this requires a non-zero 
angle $ between the laser and the motion to be cooled. 

For along-lattice motion, the effect of interaction is 
stronger, as the longest range term of C(x) (Eq. (13 1, 
Fig. [3| acts across the dipole, which with our assump- 
tions includes along the laser. At lattice spacing a = 0.5A 
(Fig. gc)), there is a high peak at ■& = (4.76 x 10" 3 
phonons per atom at d = 0, A = 10.02T), which slowly 
gets higher as N increases. This can be avoided by using 
a non-zero with a minimum of 2.05 x 10 -3 per atom 
at i? = ±0.56, A = 9.90r. While worse than the non- 
interacting case, this is much better than not optimizing. 

If the lattice spacing is allowed to vary (which can 
be done by using lattice lasers crossed at an angle in- 
stead of head-on), it again becomes possible to reach 
phonon numbers below the non-interacting minimum. 
The minimum of 1.14 x 10~ 3 per atom is reached at 
a = 0.744A,i? = ±0.05, A = 10.04r (Fig. ^d)). There 
are also local minima at a = 1.248 A and a = 1.749 A 
(Fig. [8]) , which are less deep but still go below the non- 
interacting minimum. 

In all three interacting cases, the phonon number in- 
creases a few times faster than in the non-interacting case 
as we move away from the optimum z9. However, ±5° is 



still sufficient accuracy for better than non-interacting 
cooling. 

The small waves in the contours are A^-dependent, 
decreasing in both amplitude and wavelength as N in- 
creases. We conjecture that they are an artifact of the 
periodic boundary conditions. 

We do not consider lattice spacings less than half a 
wavelength, as the Brillouin zone then contains wavevec- 
tors |p| > feoj at which Rec(p) is very small (zero in the 
N — > oo limit). Such slowly decaying excitations are ad- 
ditional Cq ~ states, so our perturbative expansion is 
not valid in this regime. 



C. Interacting case: analysis 

To study how this enhanced cooling happens, we now 
break down the phonon number into simpler parts, which 
we plot in Fig. [9j The first step in this process is to break 
the wavevector-averaged (m) ss into the individual modes 
(m pa ) ss , so all these plots are against wavevector. The 
rows (a)-(d) are the different parts, while the columns 
(A)-(E) are the 5 settings that we consider here, given at 
the top of their column. 

Fig. [9^ a) plots the single mode steady state phonon 
number 



^pa/ss 



, heating rate i?ha(p) and cooling rate 
Rca (p) 7 each of which is scaled so it would be 1 with no 
dipole-dipole interaction and the same i9, A. We see that 
all of these are strongly wavevector dependent, with all 
5 settings having both regions that are better than the 
non-interacting value and regions that are worse. The un- 
optimized cases (A),(D) have strong heating coinciding 
in wavevector space with weak cooling, while the opti- 
mized cases (B),(C),(E) mostly pair strong heating with 
strong cooling, which reduces the total phonon number 
(if H > h and C > c, then H/c + h/C > H/C + h/c). 

Fig. [9]jb) separates the heating rate into the 6 individ- 
ual diagrams in Fig. [5j while Fig. |9^c) does the same for 
the cooling rate. The dominant cooling process in all 5 
cases is phonon absorption on laser excitation (blue line 
in Fig. |9jc)), i.e. ordinary sideband cooling, as this is 
resonantly enhanced. The strongest heating process is 
usually phonon emission on decay (recoil heating, black 
line), but emission on laser excitation (blue) and the in- 
terference between them (dashed blue/black) are often 
also significant. In all these cases, the heating and cool- 
ing rates, and hence the time required to reach the steady 
state, are of the same order of magnitude as the corre- 
sponding non-interacting values. 

Finally, Fig. |9jd) plots the Fourier space dipole in- 
teraction coefficients. One effect of interaction is that 
the cooling process cannot be exactly on resonance at 
all wavevectors, as the atomic excitations have a non- 
trivial dispersion Imc(p) while the phonons have a sin- 
gle frequency v. In particular, the dips in cooling rate at 
p±ltL = ±fco, in all cases except d = across-lattice (D), 
occur because Imc(p±k L ) is large (divergent as N —> oo) 
there. (At a — 0.5A (A),(B),(E), these are the same 
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FIG. 9: (Color online) Various functions of the (one-dimensional) phonon wave vector p. Different rows are different functions, 
different columns are different settings, given at the top (lattice spacing a, cooling of either the along- or across-lattice motion, 
laser detuning and direction either standard (A = v, # = 0) or optimized (for minimum total phonon number)), (a) Heating 
rate Rh a (p) (thin red), cooling rate R ca {p) (thick blue) and steady state phonon number (thin black), each scaled to its 
own non-interacting value (dotted horizontal line). (b)/(c) Rates of the individual processes contributing to heating/cooling 
(cf. Fig. [5J: emission/absorption of a phonon either when an excitation is created by the laser (thin blue), when it hops by 
dipole-dipole interaction (thick orange), or when it decays (thin black), and interferences between these (dotted), (d) Fourier 
space dipole interaction parameters c(p) (thin blue), d a (p) (thick orange), e aa (p) (thin black), real parts solid, imaginary parts 
dotted. In all cases v — 10 T and a is the direction of motion being cooled (x for along- lattice or y for across-lattice). 
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point mod the Brillouin zone and there is hence one dip; 
at other a (C) there are two.) The along-lattice locally 
optimal spacings (Fig. M a = 0.744A, a = 1.248A and 
a = 1.749A are all close to zeros of ImC(X) (Fig. ^b)), 
which reduce this effect by removing the nearest neighbor 
(normally largest) term, reducing the size of these dips. 

These dips are accompanied by a sudden change in 
the decay rate Rec(p). The optimal spacings are also 
all near an odd number of quarter wavelengths, where 
±fco = ±7r/2a (not necessarily respectively) mod the 
Brillouin zone, dividing it into a low Re c(p) half and 
a high Re c(p) half. This gives them a particularly good 
matching of strong cooling to strong heating (a = 0.744A 
shown in Fig. |9jC), the other two are similar), which 
further reduces their average phonon number. 



We then expand both in powers of ry Q , 



-,(2) 



Pss = Pis + Pl s J + Pi 

C = £( )+£«+£( 2 ) 



(A.2) 



As in the previous section, the superscripts indicate the 
order of the term with respect to r\ a . This technique 
has previously been applied to single particle sideband 
cooling [31 E3 . It is a form of adiabatic elimination [4, 
132], but does not neglect the spontaneous decay, which 
is important as cooling requires a decay process to carry 
away entropy. 



Zeroth order 



VI. CONCLUSIONS 

To summarize, we have considered resolved sideband 
cooling in a deep ID optical lattice, and in particular the 
effect of the dipole-dipole interaction that exists between 
any dipole radiating particles. 

At the typical optical lattice spacing of half a wave- 
length, this interaction is not a small perturbation, but 
is of the same order of magnitude as the single particle 
spontaneous decay. We handle this by transforming to 
Fourier space and working in the Lamb-Dicke (deep lat- 
tice) limit, which makes the dominant dipole-dipole term 
into a dispersion of collective atomic excitations. 

We introduce an open system version of Feynman dia- 
grams, and use them to calculate the steady state density 
matrix to lowest order in the Lamb-Dicke parameter. We 
find that the steady state phonon number may be higher 
or lower than its non-interacting value, depending on the 
parameter settings, with a minimum of ~ 30% lower. 



Appendix: Direct perturbative expansion 

In this Appendix, we solve for the steady state density 
matrix by directly expanding the master equation in the 
Lamb-Dicke parameters r\ a , and verify that this gives the 
same result as earlier obtained using Feynman-like dia- 
grams. 

In the limit 77,^ — ^ (an infinitely deep lattice), we find 
that the atomic excitations are in a coherent state, while 
the phonons are stable free particles and hence can be 
in any state. As we shall see below, to obtain a unique 
phonon steady state, we need to expand the master equa- 
tion up to at least second order in i] a . 

In the following, we use the notation p ss for the steady 
state density matrix and C for the time evolution super- 
operator, i.e. 



C(pss) = . 



(A.l) 



Substituting for C from the master equation (111 and 
equating the zeroth order terms, we obtain 

-kH^ d P^~P^H^ d ]+^\ P ^) = 0(A.3) 



with if ( °> and given in Eqs. ( 23 ) and ( 27 >. To solve 



cond 

this equation, we define 



■Sk L 



C(kr) 



and s p ESpVp^k L , (A.4) 



and its Hermitian conjugate s p . This eliminates the f2 
T (o) ■ ■ 



terms from if 



H, 



(0) 



cond 



p 

U^(p) = ^Rec( P ) S >4. 



(A.5) 



In terms of Feynman diagrams, this transformation 
is equivalent to cutting off the one-line (zeroth order 
laser) vertices and their associated lines, and absorbing 
their value (which is constant, as the lines always have 
wavevector ki, and energy 0) into the coefficient of the 
vertex they are attached to (Fig. [To|a)). 

We now neglect the hardcore constraint of no more 
than one excitation per atom and treat the s' p excitations 
as bosons. This is valid in the limit of low excitation 
densities [26, 2T\ . or equivalently weak laser driving. The 
zeroth order master equation (Eq. ( |A.3[ )) is then simply 
a sum over atom and phonon modes, 



r(0) - V£ (0) 



V £ (0) 

/ j phn,pa 
P P,a 

with the atom and phonon operators given by 



(A.6) 



C 



(o) 

at,p 



£(0) 
phn,pa 



- \ c(p) spSpP - \ c * (P) P s I S p 
+Re c(p) s p ps P , 

-Wa ( b la b pa P - pbl a b pa ) ■ (A.7) 
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Since the operators C^ p and C p ^ n pa all act on differ- 
ent modes, they all commute with each other. We can 
hence choose a basis of simultaneous eigenstates of all 
these terms, which will also be eigenstates of the whole. 
Denoting the eigenstates of by p\ with eigenvalues 
A such that C^°\p\) = \p\, these are of the form 



(£) PAat,p 55 ^ PAphn.pa ■ 
P P," 

y A a t, P + y A p hn,pa 
P P," 



(A. 



with 



■'-at.p PAat,p — A at,p PAat.p • 



^phn, P a PAphn.pa — Aplm, P a PAphn,pa ■ (A. 9) 

The relevant atomic eigenstates p\ a t,p and their respec- 
tive eigenvalues A a t, p are 

PAat, P = |6p)(0 p | with A at , p = 0, 

PAat, P = |%>)(0 P | with A at , p = --c(p) , 

PAat, P = |0 p )(n p | with Aat, p = --c*(p) , (A.10) 
PAat, P = |lp)(l P | - |0p)(6 p | with A at , p = -Rec(p) . 

Here |n p ) = ^r(^p)"|0p) * s the state with n excitations 
in the mode p. The tildes indicate that these are in terms 
of the s p operators. In terms of the original s p , |0 P ) is 
a coherent state of amplitude viVn/|c(ki J )| for p = Icl, 
and a vacuum state for all other wavevectors p. 

Note that does not correspond to a unitary evolu- 
tion, and its eigenstates are hence not orthogonal. In par- 
ticular, if V is the projection onto the zeroth order steady 
states = 0, 7>(|l p )(l p |) = Vi\l v )(\ v \ - |6 P )(0 P | + 
|6 P )(0 P |) - |0 P )(0 P |, not zero. 

The phonon eigenstates are simply given by 

PAphn.pa = \™){n\ with A phn , pQ = iv a (n - m) . (A.ll) 

Since the eigenvalue only depends on the difference in 
phonon energy between the right and left sides of p, and 
phonons of all wavevectors p have the same zeroth or- 
der energy v a , these phonon eigenstates are highly de- 
generate. In particular, the zeroth order steady states 
p ss = p\ =0 are 



Pis } = < h >|0)(0|, 



„(o) 



(A.12) 



where the atomic part |0) = |0 P ) is the vacuum state of 

p 

all the s p , while the phonon part p p °^ n can be any phonon 
state which is diagonal in the total phonon energy (i.e. it 
can be a mixture, but not a superposition, of energies), 

^2 b la b Pa Pphn = hVa ^Phn ^p^P" • ( A - 13 ) 

p.Q p,a 

We see below that this degeneracy in the steady state is 
broken by the higher order terms, giving a unique p p °^ n 
when terms up to second order in r\ a are considered. 



First order 



The first order terms in rj a in Eq. ( A.l ) arc 



£(°)(pW) + £W(^)) = 0. 



(A.14) 



The unknown in this equation is pis , and solving for it 
gives 

-l 



P 



(1) _ 



£(°) (£«(/i 0) ))- 



(A.15) 



Here is formed in the same way as £(°) (inEq. (KM), 



but using the first order terms i^and an< ^ ^^(p)' a ^ so 
given in Eqs. (23) and (27). The relevant terms, i.e. 



those which ar e non -zero when applied to pis , are ex- 



tracted in Eq. (A. 21). Because has zero eigenvalues, 



the inverse is not defined on every density ma- 

trix. However, the = states pii^ are diagonal in the 



total phonon energy (Eq. (A. 13)). Every term in ei- 
ther creates or destroys exactly one phonon. These terms 
hence transform £w = states to only £( ' ^ states, 
on which the inverse is defined. This means all 

the degenerate zeroth order steady states give valid 
first order steady states. Hence, the first order terms do 
not break the degeneracy (non- uniqueness) of p p °^ n - 



Second order 



We therefore consider the second order terms in 



Eq. (A.l), which we will find do break the degeneracy 



to give a unique steady state solution. These terms are 
given by 

o = £(°)(pL 2) ) + £ (1) (pL 1) ) + /: (2) (^ ) )- (A.16) 



Substituting pis^ from Eq. (A.15) and solving for pi? 
find that 



we 



n (2) _ 



£(0) 



£(2)_£(1) U0)1 £(D (p(0)) 



(A.17) 

The term comes from processes with a single second 
order phonon vertex (the first diagram in Fig. [5]) , while 

£(!) [£(°)J £(!) comes from processes with two first or- 
der vertices. Eq. (A.17) only has a solution if the inverse 
[£(°)] _1 exists on the state it acts on, or equivalently if 
this state (in the large round brackets) has no component 
in the = subspace. This condition can be written 



= V ( J £(2) _ £(1) [£(0)1 £ (1) (p (0)) 



(A.18) 



where V is the (non-orthogonal) projection operator 
which selects al l the d ensity matrices that are of the form 



of pis' in Eq. (|A12|. Otherwise, the operator [C 



(0)1 



cannot be applied to this expression. 
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4. The V£ (2) term in Eq. (Oil) 



In the following we summarise the terms which con- 
tribute to the evaluation of the first term in Eq. (A. 18 1 



Only terms which are different from zero when applied to 
the vacuum state |0) of the atoms have to be taken into 
account, since is applied to this state. Furthermore, 
the final state of the atoms should be |0), since the pro- 



jector V selects terms which fulfill Eqs. (A. 12). We also 



only need to consider phonon operators with an equal 
number of phonon creation and annihilation operators, 
as only these keep the total phonon energy diagonal, as 
described by Eq. (A. 13). Using Eqs. (23) and (27), we see 



that the relevant terms for the calculation of the V CS 2 ^ 
term are given by 



H 



(2) 

cond / j 



mn 2 

c(kL) 

n 2 



Imc(k L ) 

VaVP "'hex khp - e Q/ 3(k L ) + e a p(k L + p) 



( b P* b lt3 + b -p a b -pl3) 



U(2) (P) = | c (k \\2 ^e afj {k L ){b pa bl p + bl pa b^ pf3 )p-2Ree a p(k L +p)(b pa pbl +bl pa pb^ p p) 



+Ree a p(k h ) p(b pa b p/3 + tf_ pQ b- p p) 



(A.19) 



where the dots stand for terms which do not contribute to Eq. (A.18), either because they are zero when acting on 
Pss\ or because the projection V removes them. Putting these together, 

f! 2 



VC^(p^) = £ 



|c(k L 



■Imc(k L )7? a 77^ k La k Lf3 - iIme Q/3 (k L ) + e a/3 (k L + p) 



2 SI 



HC_ H | c (k T ^|2 Ree "£ ( kL + P)( b P" Pss } b U + b -p« Pis } b -pP)i 



(b pa b P0 + bl pa b- pP )pW 



(A.20) 



These operators do not contain atomic operators, and are quadratic in phonon operators, with every product being 
diagonal in phonon wavevector p. In an isotropic lattice (y x = v y = v z ) a and f3 can be different, but in a fully 
anisotropic lattice (y x , v yi v z all different) we must have a = /3 to stay diagonal in phonon energy (Eq. (A. 13)). 



5. The VC (1) [£ (0) ] 1 C w term in Eq. (A.18) 

We next calculate V (jZ^ [C^] 1 (pis' ')J by applying these operators one at a time. For the first applicatic 
of D- 1 ', again using Eqs. (23) and (27), we find 



(i) = 

cond / j 2 



ilia kl^ a 



d a (k L ) - d a (p + k L ) 
c(k L ) 



{b pa + b. L )s 



S k La (b 0a + V 0a ) + ... , 



pcWp+k L 



2c(k L 



nil) (p) = Y,h^2 lmd M \p(.b 0a + bl a )-(jb 0a + bl x )p 



|c(k L ) 



(A.21) 



where . . . stands for terms which give zero when applied to the vacuum state |0) of the atoms, and hence can be 
discarded. As each term of TZ^^p) in Eq. (A.21) contains operators acting on either the left or the right of p, not 
both, these terms can be absorbed into the conditional Hamiltonian, giving 

d Q (k L ) ~d a (p + k L ) 



-(i) = ^hn 

cond / ^ 2 



c(k L ) 



(b pa + bi )s 



pa^p+k L 



hVNn 2 



Imd Q (k L ) - r\ a Rec(k L ) k ha ) (b 0a + &£ Q ) + 



(A.22) 



and TZ^(p) = 0. The master equation is in both cases exactly the same. 
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Each of the four operator terms in -ff^ond takes pis^ to another eigenstate of C^°); for example, bpaS^^pis has 
the eigenvalue — c(p + kjJ/2 + iv a , the first term from |1)(0| = Sp +k jO)(0| and the second term from the b pa . We 
can hence apply [Z^ -*] by simply inserting the reciprocal eigenvalues. Using the short hand notation 



£(0) 



fi-W JO)"* 
cond "ss I i 



(A.23) 



this yields 



STt-nf- i d a (k L )-d a (p + k h )\ ( 

> Ml vq a k La — — 



u pa —pa 

k L ) + 2iis a -c(p + k L ) - 2w a j "P+ kL 



sl 



y2 1 ^| V ^^I2 ( lm da ( kL ) _ ^ q Re c ( kL ) ^ lq ) ( b ° a ~ b < 



,(o) 



(A.24) 



This inverse is not unique: we could add in an arbitrary — state. However, doing so would have no effect on 
the final result, as it would be turned into an £,(°) ^ state by the second JCP-) then annihila ted by V . 

We now apply the second L^ x \ and finally apply V, i.e. keep only terms which satisfy Eqs. (A. 12) and (A. 13). This 
gives 



V 



. JL h (1) » + — ~H (1)t - - n<n (») 



+ H.c. 



(A.25) 



where 3 is defined by Eq. (A.24) and the Hermitian conjugate comes from [£(»)] 1 (/jffi^i^) = St. Using Eq. H, 



we find that the first term in this equation is given by 



-"cond" 



E 



h 2 n 2 



ftt h 



dff(k L + p) - rfff(k L ) 
c*(k L ) 



b-ppb 



ir/ a fc 



Lq 



d a (k L ) - d a (p + k L 
c(k L ) 



a,/3 



\ -c(p + k L ) + 2i^ Q -c(p + k L ) - 2i^ a 
x (lmd Q (k L ) - r] a Rec(k L ) k Lc ^j f&^&oa - b pb' 0a j | «- ; 
Analogously, remembering that 'P(|l)(l|) = |0)(0|, we find that 



^a|c(k L )| 4 



(0) 



(A.26) 



V 



EH 



(i)t 



cond 



E 



h 2 n 2 



d*(k L + P )-d*(k L ; 

c*(k L ) 

& -pa &- Pj 8 \ 



i?7a A;; 



Lq 



d a (k L ) -d a (p + k L ) 
c( k L) 



-c(p + k L ) + 2iz/ Q -c(p + k L ) - 2iv a 



2^ l^.m4 ^ Rec(k L ) A: L/ 3 



z^|c(k L 



(lmd a (k L ) - 77 Q Rec(k L ) fc La ) (fe 0Q ^ - &o« Ps? b < 



(A.27) 



Moreover, we obtain from Eq. (27), 



VK^(~) = - E ^ 

« c ( k Lj 

P,Oii8 



vq a k 



La 



-c(p + k L ) - 2\v a 



_ d a (k L ) - rf a (p + k L ) 
c(k L ) 

- Imc^(k L ) 



Im(i^(p + k L ) 



/jt 7) 



b-ppbl pa 



bpaPss^bpp 

c(p + k L ) + 2hv, 
-I n(°) 



-c(p + k L ) + 2iz^ Q -c(p + k L ) - 2i^ 



E 

a,/3 



hNn 4 



+bopb , 0a p 



^a|c(k L ) 

t JO) 



Im^(k L ) ( Imd Q (k L ) - Va Rec(k L ) k ha ) ( b 0a p^ b\ p - b\ p b 0a p$ - b\ a p^> b op 



(A.28) 
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Putting these three equations together, we finally obtain 



C*(k L 



c(k L ) 



Opals' - 6p Q pL J &L 6-p/3&_ pQ pL } - 6_ pQ pL } &-p/3 



E 

a,/3 



-c(p + k L ) + 2ij/ Q 
iiVfi 4 



-c(p + k L ) - 2i^ Q y 
flmd^ki,) - 77 /3 Rec(k L ) fc L) g^ flmd a (k L ) - r/ Q Rec(k L ) 
x (- WS? 6^ + 6 t^ QQ p (o) + 6 t Q p (o) ^ _ 6() ^t Q p (o) ) + h.c. (A.29) 



^ Q |c(k L ) 



Again, a and /3 are only independent in the isotropic case; in the fully anisotropic case we must have a — j3. 



6. The heating and cooling rates of the different phonon modes 



Substituting Eqs. (A. 20) and (A.29) into Eq. (A. 18), and collecting terms acting on the same phonon mode by 



relabelling — p — > p where necessary, we obtain for the phonon steady state 

= \ E [^(i^^-^^+^^ri^p^-^^t) 



+M ca0 (p)(b pa pg> 6 P _ 6 t bpa p (o)) + M c * Q/3 (p)(6 p/3 P W> bl a - p<$ b^p) 



(A.30) 



As this does not contain atomic operators, we can cancel the atomic part |0)(0|, which replaces p 8 ? by p p ° n . It is also 
quadratic in phonon operators and diagonal in phonon wavevector. 
For p ^ 0, the coefficients Mh a ^(p) and M ca p(p) are given by 



M CQ(8 (p) 



f] 2 



f "70 k L /3 



^(k L -p)-^(k L )^. d a (k L )-d Q (k L -p) 

l?7a "^Lq 



c*(k L ) 



c(k L - p) + 2ii/ a V 
fi 2 

-— — — i^[-iImc(k L ) r] a rj/3 k^k-Lp - 2iIme Q(3 (k L ) + 2e Q/3 (k L - p)], 



c(k L ) 



n 2 



. , d (k L + P )-di(k L )\ / d a (k L )-d tt (k L + P ) 

"70 «L0 17 la kha 



c*(k L ) 



c(k L + p) - 2ij/ a 
fi 2 

-— — — j2 [-iImc(k L ) r\ a r\p k^k^p - 2iIme Q(3 (k L ) + 2e Q(3 (k L + p)] 



c(k L ) 



(A.31) 



For p = 0, there appear to be additional terms 
M ha0 (O) 



rj a k La tt 2 ( f 2Imd 1 g(k L ) 

VP kh/3 - 



c(k L ) + 2h/ Q V' P P c*(k L ) 



o 2 



|c(k L ) 



-iImc(k L )77 Q 77 /3 A:LafcL / 3 + 2Re e Q)3 (k L )] 



AW(0) 



^a|c(k L ) 



(l-mdpCk L ) - rjpRec(k L ) /c L(3 ) (lmd a (k L ) - ry Q Rec(k L ) k La 



VP kh/3 



2Imrf ) g(k L )\ tt 



^a|c(k L )|4 



-i ±1110(^)^77,3 k La k U j + 2Ree a) g(k L )] 
^Im^(k L ) - ?7 /3 Rcc(k L ) A; L(3 ^Imd Q (k L ) - r7 Q Rcc(k L ) fc Lct ^) . (A. 32) 



c(k L )-2w a \- ip ^ p c*(k L ) y |c(k L )| 2 

mn 4 



However, these extra terms (the last line of each equation) are pure imaginary, symmetric in a, and of opposit e sign s 
in Mhap and M ca p. As b pa b p p — b p ^b pa = 5 a p is a number so commutes with p, such terms cancel out in Eq. (A.30), 
so Eq. (A.31 ) can also be used for p = 0. 
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(C) 
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i 

_ / * 
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+ 
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V = 



+ 



+ ooi 4 ) 



FIG. 10: (Color online) (a) Feynman diagram transformations corresponding to s p — > s p : the laser vertices are replaced by 
the vertices introduced here, the one-line (zeroth order laser) vertex being removed completely, while the hopping and decay 
vertices are unchanged. All except the one at the top right have both 1-phonon and 2-phonon versions, with the latter having 
both phonon lines attached to the same vertex. The reset vertices (the 2 in the top right, with p lines through them) always 
have at least one line on each side of p, and hence the phonon-only reset vertex (top right) has no 1-phonon version, (b) 
Diagrams for the phonon self-energy £ a (p). (c) Diagrams for the vacuum bubble V . The diagrams for £ a (p) are all on the 
left of the p line (the same diagrams on the right give E* (p)), while those for V are summed over both sides of p. As in Fig.H] 
solid orange lines are atomic excitations, dotted green lines are phonons (with either direction allowed if there is no arrow), 
blue circles are laser driving, blue stars are spontaneous decay. 



In the fully anisotropic case (three different phonon frequencies v x , v y , v z ), Eq. (A. 301 simplifies to 



= E V ^WP) - l M La(p)P^ b P abU + ReM haa (pH a pi° s ] V* 



p.a 



-M caa (p)bl a b pa - -M* aa (p)p$ bl a b pa + ReM caa (p)b pa p(°> bl 



(A.33) 



In this case each phonon mode pa evolves independently, with a heating (phonon creation) rate Re Mh aa (p) , a cooling 
(phonon decay) rate Re M caa (p) , and an energy dispersion Im M haa (p) + Im M caa (p) . 



7. Consistency with diagrammatic method 



Comparing Eq. (A. 31 1 with Eq. (29), we see that both methods give the same heating and cooling rates, i.e. 
Re Mh QQ = Rha , Re M caa = R ca . 
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The self-energy correction to the phonon propagation line, as given by the diagrams in Fig. 10 Jo), is 
S„(p) = 



O 2 

c(k L + p) - 2w a 



Vak 



La + 1- 



.d Q (k L )-d Q (k L + p) 



c(k L ) 



1 - . d a (k L ) - d a (k L + p) Imd a (k L ) 
2^ feLa 1 2c(k L ) + c*(k L ) 



iRec(kL) 



+ (same with p — > — p, ^„ — > — + 



r> 2 

c(kL) 



+ [2e QQ (k L ) - e Qa (k L - p) - e QQ (k L + p) 
~(M haa (p) + M caQ (p)) + 0(?7 4 ) . 



1 



rf a (k L ) - ri a (k L + p) 
|c(k L )| 2 
Rec(k L )\ 2Ree aa (k L ) 



c*(k L 
2Re c(k L 



c(k L ) |c(k L ) 



c*(k L ) 



The vacuum bubble (diagram with no external lines, Fig. |10[ c)) is 

n 2 



V = 2Re 



(c(k L -p) + 2LE)(ii/ a -iE) 



?y Q fcLa + i 



.d Q (k L ) - d Q (k L - p) 



,_ 1 f rf a (k L )-rf a (k L ~p) Imrf a (k L ) 
x( 2 ^fc LQ i + ^ + iRcc(k L 



c(k L ) 

d a (k L ) - d a (k L - p) 



|c(k L )| 2 



+2Re 



n 4 



i^|c(k L )| 4 
Re^M hQa (p) + 0(r, 4 ) 



i?7afcLaRe c(kij — iImd a (kL 



0(r7 4 ) 



(A.34) 



(A.35) 



where the 1/27T in the loop energy integral is a Fourier transform normalization factor, and the 2Re comes from 
including diagrams on both sides of the p line. 

Substituting these into Eq. (301 and setting p = 0, we obtain Eq. (A.13) at zeroth order and Eq. (A.33) at second 
order, i.e. the two methods agree. 
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